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Abstract 

We study sound propagation in Bose-condensed gases in a highly-elongated harmonic trap at fi- 
nite temperatures. This problem is studied within the framework of Zaremba-Nikuni-Griffin (ZNG) 
formalism, which consistent of a generalized Gross-Pitaevskii (GP) equation for the condensate and 
the kinetic equation for a thermal cloud. We extend the ZNG formalism to deal with a highly- 
anisotropic trap potential, and use it to simulate sound propagation using the trap parameters 
corresponding to the experiment on sound pulse propagation at finite temperature. We focus on 
the high-density two-fluid hydrodynamic regime, and explore the possibility of observing first and 
second sound pulse propagation. The results of numerical simulation are compared with an ana- 
lytical results derived from linearized ZNG hydrodynamic equations. We show that the second 
sound mode makes a dominant contribution to condensate motion in relatively high temperature, 
while the first sound mode makes an appreciable contribution. 
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I. INTRODUCTION 



One of the phenomena attracting attention is the superfluid dynamics in ultracold atomic 
gases. Recently, there has been renewed interest in second sound mode in superfluid Bose 
and Fermi gases l|-|5|. The existence of second sound is the most dramatic effects related 
to superfluidity in superfluid Bose and Fermi gases, which are described by Landau's two- 
fluid hydrodynamics analogous to the case of liquid 4 He {(J. These equations describe the 
dynamics when collisions are sufficiently strong to produce a state of local thermodynamic 



equilibrium [7]. In this regime first and second sound modes can be distinguished. The 
occurrence of two distinct modes is caused by the presence of both superfluid component 
and normal fluid component, which are coupled to each other. The study of ultracold 
gases in collisional hydrodynamic regime has been difficult because the density and the s- 
wave scattering length are typically not large enough. In the case of superfluid Fermi gases, 
Feshbach resonances allow ones to achieve conditions where the Landau two-fluid description 
is correct. Recent experiments have observed sound propagation in trapped superfluid Fermi 



gases with a Feshbach resonance 



We note that the occurrence of two sound modes is caused by the coupled motion of the 
superfluid component and normal fluid component. First sound is essentially an in-phase 
oscillation of superfluid and normal fluid components, while second sound is an oscillation of 
two components. This is a general feature of first and second sound, which is valid both for 
a dilute Bose gas and for superfluid 4 He. However, the detailed characteristics and behaviors 
of two sound modes are quite different in a Bose gas and in superfluid 4 He. In superfluid 4 He, 
first sound is essentially a pressure wave, while second sound is essentially a temperature 
wave. In this case, second sound is completely uncoupled to the de nsity fluctuations. The 
situation is quite different in a dilute Bose gas, as discussed in Ref. [121 ] . At very low tem- 
perature, the first sound mode is essentially the condensate collective mode and the second 
sound mode is the collective mode of quasiparticle excitations. With increasing tempera- 
ture, hybridization two modes occurs, and the nature of the sound oscillations changes. At 
higher temperature, the first sound mostly involves the noncondensate oscillation, while the 
second sound mostly involves the condensate oscillation. In the case of a Bose gas, both first 
and second sound modes are coupled to density fluctuations. Therefore, the second sound 
pole makes a significant pole to the dynamic density response function. This means that, in 



contrast to the case of superfTuid 4 He, one can probe second sound in a dilute Bose gas by 
density perturbation. 

Experimentally, sound wave in a highly-elongated trapped gas can be excited by a sudden 
modification of a trapping potential using the focused laser beam. The resulting density 
perturbations propagate with a speed of sound. This type of sound pu.se expert was 
first carried out by MIT group for a Bose gas to probe Bogoliubov sound [10]. Observed 



sound ve.oc.ty was in good agreement with theoretica. predictions 11]. In the case of trapped 
Fermi gases, first sound has been observed by the sound pulse propagation experiment [3J. 
Theoretically, sound pulse propagation in a trapped Bose gas have been studied for T = 



using Gross-Pitaevskh equation and for a normal phase using the hydrodynamic 

equation [14|. The sound propagation was also studied theoretically for a normal Fermi gas 
using the kinetic equation [15[. The sound pulse propagation in superfluid Fermi gases in 



the two-fluid hydrodynamic regime was studied in Ref. [16|. In this regime, it was shown 
that two types of sound pulses, corresponding to first and second sound, propagate with 
their sound velocities. 

More recently, sound propagation in Bose-condensed gases has been observed in Ref. 4| 
when the thermal cloud is in the hydrodynamic regime and the system is therefore described 
by the two-fluid model by using highly-elongated (cigar-shaped) traps. This experimental 
work reported evidence for a second sound mode in superfluid Bose gases, but first sound 
mode was not clearly identified. For completeness of the two-fluid hydrodynamics, it will be 
important to observe both first and second sound. 

In this paper, we study sound pulse propagation in Bose-condensed gases in a highly- 
elongated harmonic trap at finite temperatures. In order to simulate the coupled motion 
of the condensate and noncondensate components in a fully consistent manner, we use the 
formalism developed by Zaremba, Nikuni, and Griffin (ZNG) jf], [l^], that consists of a 
generalized Gross-Pitaevskii (GP) equation for a condensate and a kinetic equation for the 
thermal component. ZNG equations treat the excitations semiclassically within the Hartree- 
Fock (HF) approximation. Thus, the excitations dynamics with a thermal cloud of particles 
is governed by a Boltzmann equation for the phase-space distribution function. The coupled 
GP and Boltzmann equations include the transfer of atoms into and out of the condensate, 
which is taken together with mean-field coupling between the two components. 

We will present several dynamical simulations of sound propagation based on the ZNG 



formalism. The procedure involves solving simultaneously a GP equation for the condensate 
and a Boltzmann kinetic equation for the thermal cloud. The sound pulse is excited by the 



same manner as the experiment way J4]. In the case of a trapped Bose gas, one might think 
that the thermal density perturbations (first sound) is so small that one cannot distinguish 
small density perturbations from signal-to-noise in the thermal cloud. However, we will 
show that both first and second sound mode can be observed by a sudden modification of a 
trapping potential at intermediate temperatures. 

Since we are interested in the collision-dominated hydrodynamic regime, we have to 
simulate the system with a large number of thermal cloud atoms in order to achieve high 
enough density. However, numerical simulation of the ZNG equations for the system with 
a large number of thermal cloud atoms is very time consuming. In the present study, in 
order to save cost of numerical calculation, we derive quasi-lD ZNG equations by expanding 
the field operator in radial modes of the trap potential [18| . As shown in Ref. 18], even 
when the dynamics of Bose-condensed gases in a highly-elongated harmonic trap is well 
approximated by 1-dimensional (ID) GP equation, the momentum space of the thermal 
cloud must be treated as three dimensional (3D), because the thermal cloud atoms typically 

1 the radial trap 



171 ]. we develope 



have kinetic energy much larger than the typical energy associated wit 
frequency. Combining the work in Ref. [l8| with the ZNG kinetic theory 
the quasi- ID kinetic theory that include the degree of freedom in the radial direction. 

In Sec. II, we introduce quasi-lD ZNG equations appropriate for a highly-elongated 
3ose gas. The discussion closely follows the original approach given by Zaremba et al 
(], \v^. In this formalism, the condensate is described by a generalized quasi-lD GP equation 
for the Bose order parameter. It involves terms that are coupled to the noncondensate 



component. As in Refs Js] and fl7|. we restrict ourselves to finite temperatures high enough 
that noncondensate atoms can be described by a semiclassical kinetic equation for the single- 
particle distribution function. 

In Sec. Ill, we show dynamical simulations for a Bose condensed gas in a highly-elongated 
harmonic trap with parameters corresponding to the experiment. We also estimate the 
collisional relaxation rate which defines the two-fluid hydrodynamic regime. 

In Sec. IV, we discuss the first and second sound amplitude for condensate and non- 
condensate components separately using linearized ZNG hydrodynamic equations. In this 
section, we consider a uniform Bose condensed gas for simplicity. We calculate the relative 



weights of first and second sound mode using HF approximation for calculating thermody- 
namic various variables and compare those calculating by the dynamical simulation of the 
coupled ZNG equations. 



II. QUASI ID ZNG EQUATIONS OF A BOSE CONDENSED GAS IN A HIGHLY- 
ELONGATED HARMONIC TRAP 



We consider a Bose condensed gas confined highly-elongated harmonic trap potential. 
Our system is described by the following Hamiltonian : 

h 2 



H 



dr\ft(r,t) 



2M 



V 2 + Kxt(r) 



^ (r, t) + |^(r, t)ft(r, t)$(r, t)j>(r, t) , (1) 



with an anisotropic harmonic potential V ext (r) = ^-[u)±(x 2 + y 2 ) + cu z z 2 ]. In this paper, we 
consider a highly-elongated trap potential u z <C u±. As usual, we treat the interatomic 
interaction in the s— wave approximation with g = 4irh 2 a/M, where a is the s— wave scat- 
tering length and M is an atomic mass. In order to separate the radial and longitudinal 
degree of freedom, we expand the field operator in terms of the radial wavefunctionjl8| 



(2) 



where <p n (x, y) is the normalized eigenfunction of the radial part of the single-particle Hamil- 
tonian, which satisfies 



h 



M 



H u ± (x 2 + y 2 ) 

2M 1 2 v ' 



<p n {x,y) = E n (p n (x,y), 



(3) 



and ip n (z,t) satisfies the following equal time commutation relation : 

i> n (z,t), iJjI,^ \t) = 8(z - z')6 n!n >. 



(4) 



Using (j2J) and ([31) in (JTJ), we rewrite the Hamiltonian as 

h 2 d 2 



H=-y2 fdzft n (z,t) 



+ V ext (z) + E n 



4>n(z,t) 



2M dz 2 

+ E ^ / dz^(z,t)^l(z,t)Mz,t)M^t), 



nmlk 



where the renormalized coupling constant is defined by 



gnmlk = 9 dx dy(f)* n (x,y)(j)* m (x,y)(j)i(x,y)(l) k (x,y). 



(5) 



(6) 



The Heisenberg equation of motion for the quantum field operator ip n (z,t) is given by 



d - r~ 
ih—J) n (z,t) = il) n (z,t),H 
at L 

h 2 d 2 



2Mdz 2 



+ V ext (z)+E n 



$ n (z, t) + J2 9n m iA( z , t)#i(z, t)$ k (z, t). (7) 



rnlk 



In order to deal with the Bose broken symmetry, we separate out the condensate wave- 
function from the field operator as 



$ n {z,t) = ® n (z,t) +1p n (z,t), 



(8) 



where the condensate wave function is defined by Q n (z,t) = (4f n (z,t)). The equation of 
motion for <3> n can be obtained by taking statistical average of ([7]) : 



ih^<f> n (z,t) 



h 2 d z 
2Mdz 2 



+ V ext (z)+E n 



$ n (z, t) + J2 (9nkkin k c + 2g nkkl n k ) t) 

kl 



+ J29nkkim k ^* l (z,t) +J29nmkl(^L( z ^) 1 Pk(z,t)lpi(z,t)}, 



(9) 



mk 



mkl 



where = &* k (z)& k (z), n h = (Tpl(z,t\d>k{z,t)). In Eq. ([9]) we have neglected the anomalous 
average (ipk(z, t)ip k (z, t)), as in Ref. [6J. Moreover, we have assumed that the off diagonal 
terms of the noncondensate density, (V'/tVvH^ are smau an d thus can be neglected. 

This assumption is expected to be valid in the case oo z <€i uj± In addition, in the case 
where u z C wi, the contribution from higher radial modes to the condensate wavefunction 
is negligibly small Therefore, we will henceforth approximate $ a = $5 a> o. With these 
approximations the generalized GP equation ([9]) reduces to 

h 2 dz 



2M dz 2 



+ V ext (z) +E + goooon c (z, t) 



+ 2gokkon k {z, t) - iR{z, t) 
k 



$(z,t). 



(10) 



Here the source term R is given by 



R(z,t) 



12 



2n c (z,t) 

with r i2 = -21m (j2nmk9on m k®*{^Uz,t)i/j m (z } t)^ k (z } t)}y n c (z,t) = \${z,t)\ 2 and n k 
$l{z,t)$ k (z,t)). 



(11) 



We now turn to the dynamics of the noncondensate. The physical properties of interest 
are in principle defined by the following equation of motion obtained from ([7]) and (Q: 



ih—ip n {z,t) 



2M 



+ V ext (z) + E n 



ip n {z, t) + 2 9nkkin k {z, t)ipi{z, t) 

kl 

-2 9nkkin k (z, t)i>i(z, t)+J2 9nooi®(z, t)$(z, t)$i{z, t) 

kl I 

+ J2g n mko$(z,t) i! m (z,t)tjj k (z,t) -m k (z,t)5 mk 



ink 



+2$>nmfcO$0M) \ft m (z,t)4> k (z,t) -fl k (z,t)5, 



ink 



ink 



+ J2gnmkl ^(2,t)^fc(z,t)V>/(«,t) - {$ m (z,t)$k(z,t)i>l(z,t)) , (12) 



rnkl 



where n = n c 5o yk + n ■ It is convenient to define the time evolution of ip n (z,t) by 

$ n (z,t) = S 1 (t,t )ijj n (z,t)S(t,t ), 

where the unitary operator S(t,t ) evolves according to the equation of motion 

ih-^-S(t,t ) = H eS S(t,t ), 
at 

with S(to,to)=l. The effective Hamiltonian in ( Fl4|) is given by 



H e s — Hq + H , 
H = H[ + H' 2 + H' 3 + i?4, 



(13) 



(14) 



(15) 



where the various contributions are defined as 

h 2 d 2 



z,t) 



2Mdz 2 

U n (z,t) = V ext (z) + 2j29nkk n n k {z,t), 



+ E n + ft n U n (z,t))^ n (z,t), 



H[=J2 [dz (L^l + L&n), 

n J 

Li = ^2g n mki 2n fc <M (M 5 mfc + m k ^*5 0j i5 mk + (ft^^i) 



rnkl 



900kl 



kl 



dz 



H' 3 = E 9omko f dz + $*i>li> k ipi 



rnkl 



H' 4 = J dz {^Y^^ln^i - 2g 



nnkin n i)\^l 



nmkl 



(16) 
(17) 

(18) 
(19) 

(20) 

(21) 

(22) 



The expectation value of an ordinary operator defined in terms of ip n and ip^ is given by 

(6(t)> = (0) t = Trp (t )O(t) = Trp(t, t )6(t ), (23) 
where p(t,t ) = S^(t,to)p(to)S(t,to) satisfies the following equation 

H c s, p(t, t ) 



dp(t,t ) 



dt 



(24) 



Our ultimate objective is to obtain a quantum kinetic equation for the noncondensate 
atoms. We define the Wigner operator as 



fniPz, Z, t ) = j dz'e^'^l \z + -, toj *Pn lz ~ t 

The Wigner distribution function is then given by 

f n (p z , z, t) = Trp(t, t )f n (p z , z, t ). 
The equation of motion for / is obtained by using Eq. (E 



(25) 



(26) 



dUv ^ t] = i-Trp(t, to) [f n (p z , z, t ), H Q (t)) + i-Trp(t, t ) [/ B (p„ z, t ), Jf(t)]. (27) 

With the assumption that U n (z,t) varies slowly in space, we then have 

[f n (p z ,z,t ),H (t)} ~ -—p z —f n ( Pz ,z,to)+ih—U n (z,t)—f n (p z ,z,t ). (28) 

The second right hand side of the right hand side of Eq. (128]) represents the effect of collisions 
between the atoms. As we show in Appendix A, the collision integral is the sum of two 
contributions: 

dfn 



dt 



coll 



C l2 [fn]+C 22 [f n ]. 



(29) 



Thus we obtain 



df n ( Pz ,z,t) TP ± d_ f , zt) _± U n (zt) _?_ f( zt) 

= C 12 [f n ] + C 22 [f n ]. 

The C\ 2 collision integral is defined as the contribution from the H' 3 perturbation 
Ciat/n] = -iTrp(t,to) f n (p z , Mo), H' 3 (t) 



(30) 



X] 9n'm'k'0 n c X] 'H'' ' ' '' ' 2" C 3 )j "/'■•■ + /' i /' J + ;' .; 

n'm'k' Pzi,Pz2,Pz3 

(<Wp*l<W - <$p z ,p, 2 <W ~ ^p z ,p z3 ^nk') [(1 + /f ^ _ fl + )( 1 + fs )}; 



5, 



(31) 



2 

where the local HF single-particle energie is e l = ^fy + E x + U l (z,t) and e c = \x c + 



The local condensate chemical potential ji c is defined by 



h 2 



He = 







2MJn c (z,t) dz 



\Jn c (z,t) + V cxt + E + g oo n c + 2g kkon k , 



(32) 



and the condensate velocity is given by v c = jq-§^6(z,t) with &(z,t) = <Jn c (z,t)e ie ( z,t \ The 
source term R is directly related to the C 12 collision term 



R(z,t) = 



h 



E 



dp z 



C 12 [fn}. 



2n c (z, t) ^ J 2nh 

Similarly, the C 22 collision is defined as the H' A perturbation, which is obtained as 



(33) 



C 22 [/n] 



-iTrp(t,t ) \f n (p z ,z,t ),H' 4 (t) 



~k' 



~m 

C' 



' - el) 



— n E 9n'm'k'l' E 

n'm'k'V Pzi,Pz2,Pz3,Pz4 

^Pzl+Pz2,Pz3+Pz4\^Pz,Pzl^nn' + 8p z ,p z2 3nm' ~ ^p z ,Pz3^nk' ~ ^ 

[(i + /r')(i + /ft/f/f - /r'/ 2 fc '(i + ff)(i + /!')]■ 



(34) 



We refer to Appendix A for detail derivations of the collision integrals Ci 2 and C 22 . 

In summery, we have obtained a coupled set of equation of motion for the condensate 
and noncondensate as follows: 



ih^$(z,t) 



h 2 dz 



2M dz 2 
+ E 2gokkon k (z, t) - iR(z, t) 



V ext (z) + E + g 0000 n c (z, t) 



®(z,t), 



df n (Pz,z,t) . Pz d d d 

dt + ~M~dz *' t] - d~z U ^ l) dp-J n{Vz ' *' t] 

= C 12 [f n ] + C 22 [f n }. 



(35) 



(36) 



The condensate is described by a quasi-lD GP equation for $(z,t). The noncondensate is 
described by a quasi- ID kinetic equation for the distribution function f n (p z , z, t). Here n is 
the radial mode index. Different radial mode are coupled through the mean-field interaction 
as well as collisions. 

Before closing section, we give equilibrium solution of the coupled ZNG equations. The 
equilibrium solution for the condensate wavefunction is given by $ (z,i) = $ (<2)e~ lMc0 */ a , 
where $o(z) satisfies 
h 2 dz 

+ V ext (z)+ E + g ooon c0 (z) + Yl 2 9okkon k (z, 

k 



2Md 2 



%(z)= f i c0 %(z). (37) 



9 



Here //o is equilibrium chemical potential. The equilibrium distribution function is given by 
the static equilibrium Bose distribution 

f ° n{Pz ' Z) = exp{A,[p 2 /2M + t^(z)-^ o]}-l' (38) 

where 0q = is the inverse uniform temperature. The trapping potential is augmented 
by the HF mean-field U£ = V cxt (z) + 2J2k QnkkrJiQ^z). The coupled equations ( |37|) and ( 1381) 
must be solved self-consistently. 



III. DYNAMICS OF FIRST AND SECOND SOUND IN A BOSE CONDENSED 
GAS 

Using the quasi-lD ZNG equations derived in the previous section, we study sound pulse 
propagation excited by a sudden modification of a trapping potential. The numerical pro- 
cedure for calculating ZNG equations closely follows that described in Ref. jf], Q. The 
dynamics of the thermal cloud is calculated by using iV-body simulations The dynam- 
ics of the condensate is determined by numerically propagating the GP equation using a 
split-operator fast Fourier transform (FFT) method. The numerical method is described in 
detail in Appendix B. 

We take the physical parameters from the experiment of Ref. [l], which reports the 
observation of second sound propagation. In this experiment, total number of 23 Na atoms 
N = 1.7 x 10 8 , radial trap frequency u; ra d/27r = 95Hz, and the aspect ratio w ra d/w ax ~ 65. In 
this situation, one has a high density cloud of n Q ~ 10 20 cm -3 . At the lowest temperatures, 
the BEC has a radial TF radius of roughly 22 /im and an axial TF radius of 1.4 mm. The 
number of test particles is ten times the actual number of thermal atoms in order to minimize 
the effects of a discrete particle description. 

We first consider equilibrium solutions (1371) and (|38|) for these experimental parameters. 
In Fig. [TJ, we plot the condensate fraction N c /N as a function of the temperature. We see 
that the transition temperature for the Bose-Einstain condensation is given by T c ~ 350nK. 
In Fig. El we plot the equilibrium density profiles of the condensate and noncondensate at 
T = 176.5 nK(~ 0.5T C ). For comparison, we also show equilibrium density profiles obtained 
from the full-3D HF calculation, i.e. without making the quasi-lD approximation, in Fig. El 
The differences between with and without the quasi-lD approximation are only a few %. 

10 



This confirms that our quasi- ID treatment can well describe the highly-elongated system. 
We note that in order to obtain reasonable results, we must take large enough number of 
radial modes so that E n > k B T. For example, we took about 1000 radial mode for the 
calculate of Fig. [2J 
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FIG. 1. The condensate fraction N c /N as a function of the temperature. 
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FIG. 2. Axial density profiles in the equilibrium state. Gray lines are solutions within the quasi-lD 
approximation. Black line are equilibrium solutions without making the quasi-lD approximation. 

We now consider density disturbance by a sudden modification of the external potential 
generating pulse propagation. Here we set the external potential SU(z,t) = Ae~ az2 9(—t) 
with A ~ 0.5yU and a ~ 250 v /5 ~ with Zq = ^ Mu/h. A localized potential is applied at 
t < 0, while it is turned off at t = 0. This situation can be described as presence of a 
localized potential aimed at the center of the trap, which acts as a repulsive trap. Turning 
the potential suddenly off causes a local dip of the BEC density. This perturbation splits up 

11 



in two waves propagating symmetrically outward, both with half the amplitude of the initial 
perturbation. A schematic representation of the excitation procedure is shown in Fig. |3j 
The axial density profiles, shown in Fig. H]for various propagation times, clearly shows that 



FIG. 3. Schematic representation of the excitation of a sound wave, where the trapping potential, 
height and width of the perturbation are roughly on scale. 

density dips corresponding to two sound modes travel with their sound velocity. The faster 
sound pulse corresponds to the first sound, while the slower sound pulse corresponding to the 
second sound. We see that the both depth of the first and second sound dips at T ~ 0.5T C 
are sufficiently large for the experimental observation. For comparison, we plot axial density 
profiles of condensate and noncondensate separately for various propagation times in Fig. 



The axial density profiles at T = 88.2 nK(~ 0.25T C ) and T = 264.5 nK(~ 0.75T C ) are 
shown in Figs. |6] and [7J Compared to Fig. [5] with Fig. [6] and Fig. [7J we see that the first 
sound mode is dominant at low temperature, while the second sound is dominant at high 
temperature. 

We note that it is difficult to observe sound propagation in a noncondensate thermal 
cloud because the thermal density perturbations is so small that one cannot distinguish 
small density perturbations from signal-to-noise in the thermal cloud. Nevertheless, at the 
intermediate temperature T ~ 0.5T C both first and second sound pulses appear in the total 
density. In this regard, we note that the analysis of Ref. [41] was based on the assumption 
that the condensate motion is always dominated by second sound at all temperatures. How- 
ever, the calculation shows that the condensate motion is dominated by first sound at low 
temperatures. Therefore, it is possible that the experimental result of Ref. at T < 0.5T C 
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FIG. 4. Axial density profiles of condensate and noncondensate for various propagation times at 
T ~ 0.5T C . 
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FIG. 5. Axial density profiles of total density for various propagation times at T ~ 0.5T C . 

may have observed first sound. It may therefore require a careful analysis of the experi- 
mental data in the crossover temperature regime T ~ 0.5T C in order to identify two sound 
modes. 

Let us now examine validity of two-fluid hydrodynamics in our system. The existence 
of first and second sound is predicted by Landau two-fluid hydrodynamics, which is valid 

13 





FIG. 6. (a) Axial density profiles in the equilibrium state without perturbation (i. e., 5U = 0) at 
T = 88.2 nK(~ 0.25T C ). (b) Axial density profiles at propagation time t/to = 30. 
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FIG. 7. (a) Axial density profiles in the equilibrium state without perturbation (i. e., 5U = 0) at 
T = 264.5 nK(~ 0.75T C ). (b) Axial density profiles at propagation time t/to = 30 

when collisions are sufficiently strong to produce a state of local thermodynamic equilibrium 
This requirement is usually summarized as lot 1, where u is the frequency of a 
collective mode and r is the appropriate relaxation rate. In a trapped Bose gas, a relevant 
relaxation time is r 12 relaxation time associated with the C\2 collisions jf], 1^] , which describes 
equilibration between the condensate and the thermal cloud. The relation time T\2 is given 
by 

l/ri^E 11 !^, (39) 



where 



n 

rSU] = E I d P^ I dp* [ dp z ,6(e c + ef - ef - ef) 

„W nh J J J 

X8(pcz + Pz2 - Pz3 ~ Pz4)(Snn> ~ Sum* ~ 6 nk >) f%' \l + ff)(l + /|'). (40) 
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In Fig. |H1 we plot the equilibrium local collision rate 1/tu in a trap as a function of the 
z distance at T ~ 0.75T C . This collision rate again has a maximum at the edge if the 
condensate, falls off rapidly beyond this point, being proportional to the condensate density 
n c . This figure shows that Wo r i2 < 1 in the whole region of the condensate, where Uq = ui£, 
£ = ^ 2 ^ gn being the healing length. For the trap parameters given above and in the 
temperature range 0.25T C < T < 0.7T C , we found that the equilibrium local collision rate 
satisfies u; t 12 < 1 in the whole region of the condensate. Thus, the sound propagation 
experiment is well within the hydrodynamic regime [ijj]. For comparison, we also calculated 
I/T12 without making the quasi-lD approximation = J dxdy^p?/ J dxdyfio [6j. The 

differences between with and without the quasi-lD approximation are only a few %. This 
also confirm that our quasi- ID ZNG equation describe dynamics of the highly-elongated 
system quite well. 




FIG. 8. The equilibrium local collision rate I/V12 in a trap (in units of a frequency Wo = ui£ with 
the speed of first sound u\ and healing length £ = —?===) as a function of the z distance (in units 
of the harmonic oscillator length) at T ~ 0.75T C . 



In the collisionless regime, i.e. ooqTu > 1, the first and second sound cannot be excited, but 
only Bogoliubov sound can be excited by a sudden modification of a trapping potential. In 
Fig. Owe show the sound pulse propagation at total number of atoms iV = 10 5 corresponding 
to the density n = 10 14 cm" 3 at T ~ 0.5T C . We see that the dynamics of the cloud is quite 
different from Fig. |5j Clearly, only the Bogoliubov sound propagates. This situation is 
similar to the case of MIT experiment in Ref . [uj] . 
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FIG. 9. Axial density profiles of condensate and noncondensate for various propagation times for 
iV = 10 5 at T ~ 0.5T C . In this case, the system is in the collisionless regime. 

IV. COMPARISON WITH LINEAR RESPONSE SOLUTION OF ZNG HYDRO- 
DYNAMIC EQUATIONS FOR A UNIFORM GAS 

In this section, we derive an analytical expression for the amplitude of first and second 
sound pulses for a uniform gas. We limit ourself to the hydrodynamic regime. General 
expressions for pulse amplitudes have been derived using Landau two-fluid hydrodynamic 



equations in Ref. 



16]. Here we instead use ZNG hydrodynamic equations, which allows for 



direct comparison with the simulation results. We start with the linearized ZNG hydrody- 



namic equations for a uniform gas 20| : 

d5n 



w -n (V-Sv n ) + 5T l2 , (41) 
at 

or 

Mn -^ = -V5P - 2gn V(5n + 5n c ) - h V5U, (42) 
85 P 5 ~ 2 

— = --P (V ■ 8v n ) + -(/i c0 - U Q )6T 12 , (43) 

^ = -n c0 (V-5v c )-5r 12 , (44) 
M^ = -V(6fr-6U), (45) 

where 

5 fi c = gSn c + 2g5h . (46) 
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The expression for SF 12 is given by 



5T 12 [f] = - ^ 0?7,c0 ^diff, /idiff = /2 - /i c , (47) 
T12 



HcO = gn C Q + 2gn, U = 2g(n c0 + n). (48) 

In the above equations, we have explicitly included the time-dependent external perturbation 
5U(r,t). To solve the linearized hydrodynamic equations, we introduce velocity potentials 
according to 5v c = V0 C and 5v n = V0 n . In terms of these new variables, the equations for 
the condensate and the equations for the noncondensate can be combined to give 

c^2 1 (J c) 

M ~&F = ^ coV2 ^ c + 2 9n Q V 2 <p n + ^-5/idifr - -^SU, (49) 

= ( ?5 + 2^1 V 2 n + 2^ cO V 2 c - ^^<W - (50) 

Here <5r 12 has been expressed in terms of Sfi^is using ( 14T|) . The equation of motion for <5/Xdiff 
is given by 

a. = o^coV 0„ - #ri c0 V (p c , (51) 

ot 3 r M 

where the relation time r M associated with the chemical potential difference is defined via 

1 _ Pogn c o ( \h + 2gn n c0 + \^gn\ \ _ f3 gn 
T n ~ T i2 \ |7o^o + \gn\ ) ~ r 12 cr H ' 

As discussed in Refs. [if]], the liner response to the pulse perturbation can be described 
in terms of the density response function x„ n (q, w). We will thus calculate x„ n (q, ui) by 
considering the external perturbation that excites plane wave SU = 6U qju e i(n ' r ~ ut) ■ Therefore 
we look for the plane-wave solutions Ci „(r, t) = Cj „ i q iW e l( ' q ' r ~ tJ '\ In this case, (ISTl) reduces 
to 

5/idiff = ~ — — gn c0 ((pc ~ 770n) q 2 - (53) 

1 — iujt^ \ 3 / 

Substituting this result into ( 1491 and (15T?|) . we are left with two coupled equations for the 

superfluid and normal fluid velocity potentials: 

2 / an \ 2 
Mu (p c ^ u = gn c0 1 - ; q <p c ^ 



1 - iuJT fx/ 

,-. . (Th n c0 

+2gn 1 + — ; - — 

3(1 - iUT^) n 
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q ^n,q, w - iuSU^, (54) 



and 



Mu 2 



5Po 
3n n 



2^c 



2a 



H 



n 



cO 



9(1 - iwr M ) n 2 . 



+ 2c/n 



cO 



1 + 



OH 



3(1 - zwtJ n 



Taking the limit UTn — > of these coupled equations, we obtain 



y ri,q,o; 



#n c0 (l - o-£r)g c , q ,o, + 2fifn (l + ^ j g n>q 



iu5U. 



5P o , 2°Hn 2 c0 



— + 2gn 1 
3n \ 



9^ 



y ?i,q,a; 



+ 2#n c0 1 + 



3n / 



<7 



y c,q,w 



iujbU. 



q,w 



It is useful to rewrite (|56|) and (1571) in a simple matrix form as 

/ u 2 -v 2 q 2 -v 2 2l q 2 \ I . 

I = —iujbU, 

\ -vf 2 q 2 u 2 -v 2 q 2 ) y<f> 
where we have introduced new velocities 

.2 



'12 



M 

2gn c0 
M 



1 + 



3n / 



3n / 
^- + 



5P , 2^0 A 2<?Hn 2 c o 



9n 2 n 



1 3Mh M \ 

We note that these new velocities are related to the first and second sound velocities 
u 2 through 



2,2 2,2 22 22 22 

u 1 +u 2 = v 1 +v 2 , u x u 2 = v x v 2 - v X2 v 21 . 



Solving fl58l) . we obtain 

/ j. 



\ 



'n,q,u; 



-iujSU, 



(u 2 — u 2 q 2 )(u 2 — u 2 q 2 ) 



uj 2 — v\q 2 + v 21 q 2 
\ v i2<f + - v 2 q 2 



Using ( 1531) and Taking the limit wr^ — > 0, (H5l) and ( HT1) reduce to 





- 92 f 







rz c0 (l + <r H ) 



t 



\ -v H n c0 n (l + |^; 
Using the solution flBTj) in the expression fl62l) . we obtain 



^c,q,u 



5n r 



q 2 6U, 



q,cj 



^ M(co 2 -u 2 q 2 )(u 2 -u 2 q 2 ) 



( 



x 



n c0 (l + \cr H )u 2 + {n c0 (l + a^)(-^ 2 + T&) + |^n c0 K - ^ 2 2 )}? 2 
^ n (l - V 2 + Ko^i - t&) + no(l + + v 2 2 )}q 2 
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The solution can always be written in terms of the density response function, defined as 




SU. 



Xr 



\ Xf 



(64) 



with 



Xr 



Xf 



q 2 n c0 (1 + \a H )u 2 + {(1 + <Jh){-vI + vli) + f ggfol ~ 

M {u 2 - u\q 2 ){u 2 - u 2 2 q 2 ) 

q 2 n (1 - + {tf^K - v 2 21 ) + (1 + |^)(-^ 2 + ^ 2 )}g 2 



(65) 
(66) 



M (to 2 - u\q 2 )(u 2 - u\q 2 ) 

In the case of the sound propagation experiment, a localized potential is applied at t > 0, 
while it is turned off at t = 0. This situation can be described as SU(r, t) = 5U(z)6(—t) 
[if]]. In this case, the density fluctuations at t > is given by 



Sn c (z,t) : 
5h(z, t) 



2tt 2 
I 



dco5U(q) , Xncn . e ig2 ~^ (*>0), 
v 7 (u? + Z7?) v ; 



/ du5U(q)- 



Xf, 



jiqz—iujt 



{t > o), 



(67) 
(68) 



2tt 2 J J (w + in) 

where xj^jq.w) = Im X« c n(q, w + 177) and xL(q^) = ImXn«(q, w + 177). From (JB2J and 
(!68|) . we obtain 

5n c (z, t) = W?° [5U(z - uit) + 5U(z + Ul t)\ + [5U(z - u 2 t) + 5U(z + u 2 t)} , (69) 
5n{z, t) = W? [5U(z - u x t) + 8U(z + Ul t)\ + W? [5U{z - u 2 t) + 5U(z + u 2 t)} , (70) 



where the amplitudes of the sound pulse are given by 



n c0 (1 + \a H )u\ + (1 + ch){-v \ + v^) + \a H {i 



J 12) 



2Mu{ 

n c0 (1 + \a H )u 2 2 + (1 + cr H )(-v 2 + v 2 21 ) + \a H {v\ 



J 12) 



2Mul 



n 



no 



J 2\. 



+ (1 + 



2 n c0 

3 n 



)(-V 2 2+V 2 12 ) 



2Mu{ 
n 



"7 



1 - ¥h^)u\ + ^a H (v 2 - v 2 21 ) + (1 + m{-v 2 



no 



3 ftp 



2Mul 



(71) 
(72) 
(73) 
(74) 



We estimate the interaction parameter for a uniform gas corresponding to the experiment 



ofRef. 



4] from the average density of the trapped gas, and obtain 



na 



0.07 and thus 



k B T° 



0.5, were T c ° is the BEC transition temperature of an ideal Bose gas. In Fig. [TTJl we plot the 
first and second sound velocities as a function of temperature within the HF approximation, 
and compare with the sound velocities deduced from ZNG simulations discussed in the 
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previous section. We emphasize that both first and second sound velocities obtained from 
ZNG simulation show good agreement with HF approximation. This confirms that our ZNG 
simulation well describe the two-fluid hydrodynamics. 

1.2 
1 

0.8 

o 
O 

> 0.6 
0.4 
0.2 


0.2 0.4 0.6 0.8 1 

T/Tc 

FIG. 10. First and second sound velocities as a function of temperature. 

We now compare the pulse amplitude obtained from simulation results in the previous 
section with the ZNG hydrodynamic results for a uniform Bose gas. The amplitudes of 
first and second sound W™ c and are obtained by taking average of subtracting the 
unperturbed density profile from perturbed ones. In Fig. [Til the first an d second sound am- 
plitudes for condensate and noncondensate components obtained by the simulation of ZNG 
equations (Eqs. ffTUj) and f l30|) ) and the results f J7Tj) -f TMl) which calculated by the linearized 
ZNG hydrodynamic equations. We see that the simulation results are consistent with the 
analytical results. 



ZNG simulation 
- HF approximation 



First sound / 




V. CONCLUSION 



In this paper, we have discussed sound propagation in Bose-condensed gases in a highly- 
elongated harmonic trap. In order to consider the situation of a highly-elongated harmonic 
trap, we derive quasi-lD ZNG equations. Using these equation, we show the several dynam- 
ical simulation with the same parameter as experiment on second sound. We showed that 
both first and second sound mode can be observed by a sudden modification of a trapping 
potential at intermediate temperatures. We also found that the thermal density perturba- 
tions is so small that one cannot distinguish small density perturbations from signal-to-noise 
in the thermal cloud. 
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FIG. 11. The first sound amplitude Wi = W? c + W? and second sound amplitude W 2 = W£ c + W% 
as a function of temperature, where W = Wi + W 2 = W™ c + W£ + W 2 " c + W.f . Lines show the 
results from the self-consistent HF approximation. Data points are results from the simulation 
solving ZNG equations. 



We also derived expression for the pulse amplitude of condensate and noncondensate 
components in a uniform Bose gases using linearized ZNG hydrodynamic equations. The first 
and second sound amplitude obtained by dynamical simulation are consistent with the results 
calculated by the linearized ZNG hydrodynamic equations. This also confirm that the system 
we considered in this paper is well described by the two-fluid hydrodynamics. The quasi-lD 
ZNG formalism developed in this paper is very useful in analyzing the finite temperature 
dynamics of highly-elongated Bose-condensed gases. In a separate paper, we will study the 
collective modes of a highly-elongated Bose gas at finite temperatures. We hope to stimulate 
further detailed experimental examination on the dynamics of Bose condensed gases at finite 
temperatures. 
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Appendix A: DERIVATION OF COLLISION INTEGRALS 



In this Appendix, we give a detailed derivation of the expressions for the collision integrals 
given by (15T1) and We closely follow the approach of Refs. jf], \v\ . 

According to the time dependent perturbation theory, the expectation value of an arbi- 
trary operator 0(t) made up of some combination of non-condensate field operators can be 
expressed to first order in H' as 

(d) t = Tr{s (t,t )d(t )So(t,t ) 

rft / 4 t (t / ,to)[^(t,t / )0(to)5o(t,t , ),^ , (^)]5o(t / ,to)). (Al) 



ft Jt 

The three-field correlation function is given by 

= -~Trp(* ) / < ^(t\t o )[4 t (t,t0^(^^o)^ m (^,to)^(^^o) 
n Jt 

xS (t,t'),H[(t') + H! 3 (t')]S (t',t ) 
= (^t( M )^ m (^)^ t ))(i) + (i,}Xz,t)^ m (z,t)Mz,t))^, (A2) 

where (• • •) ( '** ) denote the contribution from H[. In this case, contributions from the other 
terms in H' can be shown to vanish. The evaluation of correlation functions is facilitated by 
two key assumptions: The effect of H[(t') in the interval to < t' < t is essentially a collision 
process, which occurs on a time scale much shorter than all other time scales in the problem, 
and the hydrodynamic variables vary slowly in space and time. It is sufficient to use 

n n c {z',t) ~ r%(z,t) ~ n c (z,t), h n (z',t) ~ h n (z,t), 

6(z\ t) ~ 6(z, t) - 6 -^(t' -t) + y _ z)i (A3) 

and 

S Q {t,t') ~ e-^oWft-f), ( A4 ) 
Introducing the Fourier transform of the non-condensate field operators according to 

4>n(z,t ) = J2 a Pz,ne ipzZ - (A5) 
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The contribution from H[ to the commutator is given by 
[Sl(t, t')4>l(z, t )tp m (z, t )jn(z, to)So(t, t'), H[{t% 

Tl'fc'J' PzlPz2Pz3Pz4 

n'V PzlPz2Pz3Pz4 

x ^Pz4Pzc[ a l zl n a Pz2rnSp z3Pz4 5i/ > i + aj, zin ap z3 i(5p z2 p z4 (5i/ im ], (A6) 
where we have defined the condensate momentum p zc = mt> c . We thus obtain 

<$(z, t)j> m (z, t)Mz, *)>(i) = -^Trp(to) £ df'Sj^, *o) 0$,(*, *otym(*, *o)^i(^, ^o) 

x S (M'),#i(0]So(i',to) 

= 2iE5nwn ,, '(M)^ M 

n'J' 

E p -«(Pzc+Pzl-p Z 2-Pz3)z/?»X 
C "Pz4Pzc 



PzlPz2Pz3Pz4 



X 



•'tn 



+ { a p zl n a Pz3l)3pz2Pz4. S v , m ], (A7) 



where 



(aLn<w)f = Trp^o)^^,^)^^^^^^,^) ^ ^ 1 ~' €iW ~^\al tin a Ptim ) to . (A8) 

We now assume that the initial statistical density matrix p(to) has the form appropriate for 
the HF Hamiltonian 

( a p z in a Pzim)t() — 3pzl,Pz2^nmf {Pzli z it)i (-^-9) 

where we have assume that off-diagonal contribution n ^ m is neglecting small. In fact, in 
equilibrium, f nm (n ^ m) < 10~ 5 and J2 n ^m fnm < 10~ 3 . Similarly, we can obtain 

(fy(z,t)j; m (z,t)ii>i(z,t))(3) 

= ~Trp(t ) / t d* , 5j(*',*o)[5j(t,* , )^(^to)^m(^*o)^(^,*o) * 5 (i, f), H' 3 (t')]S (t' , t ) 
n J to 

= -2 E gn'm'VonWe" E e^ 2+ ^™-^> 

n'm'fc' Pzl,Pz2Pz3 



z3 
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X [8p2,p' 2 8p 3 ,p'Jmm'fak'{al zl n a p' zl ri}t' + &p^p'j>mm! ( a Lm a p' zl n')t' ( a l z3 k a p' z3 k')f 
+Sp 3 ,p' 3 6ll' ( a lzW a p' zl n')ti ( a l z2 m a p' z2 m')t' ~ Sp^Snn' ( a l z2 m a p' z2 m')t' ( a Pz3 k a p' z3 k') V 
+^p 2 ,p' 2 ^mm'(al zl n a p z 3k)t'(al' z3k 'ap' zl n')t' + Spi,p{Sn,n'((4, z in a Pz2m)t'(ap' z2m 'ap' zl n')t' '• 



The last two terms reduce to 



(A10) 



v' zV v' z2 v' z3 



E ^2,pi^+pi 1 ,p' z2 +P^(^mm'(ap zl n^3&}t'( a p;3ifc' a p' J! in , }t' 

+ ^pi^n,n'(al zl n a p z2 m)t'(ol^ 2m/ a p ^ in ')t') 
E dp zc ,p' z2 (dpz2,p' z2 dm,m' ( a Pzl n a Pzik}t + ^p z3 ,p' z2 ^k,m' ( a p zl n a Pz2m) t 



(All) 



Pz2> 



We now see that the last two terms in this equation exactly cancel the contribution from 
H[. We thus obtain 

(ip j n (z,t)^ m (z,t)^ k (z,t)) t = -i2vrs( nmfc onyV e 



E 



Pzl,Pz2,Pz3 L 



S(e c + e? - e? - e£) + -P- r- 



5pz C+ pzupz2 +P z 3 if?a + / 2 m )(i + m - a + m m m (A12) 

with fl = f l (p Z i, Zi, t). In addition, we have treated the system as locally homogeneous, with 
the consequence that the local HF single-particle energies, e l = + E\ + U l (z,t). Using 
( 1A12I) . we obtain the Cn collision integral 



Ci 2 [f n ] = -iTrp(t, t ) f n (p z ,z,t ),H' 3 (t) 



E 9r, 



In'm'k'Onl^e 18 ^ E elqzZS (Pcz+Pzl,Pz2+Pz3)e 
n'm'k' <?z Pzl,Pz2,Pz3 

J „ .„ ..\._X -X ./J 



ipzcZ 



[^Pzi,Pz+qz/2^mn' ( a Pz -q z /2n a Pz2m'0'p z3 k'}t &p z2 ,p z -q z /2&nm' { a Pzl n a Pz-qz /2m«p z3 fc') t 
~^Pz3,Pz-qz/2^nk'{a\ )zin a'p z2 m'0'p z -q z /2m)t ~ h.C.) 

^ E glm'k^c E [ s ( e c + e?' - 3P' - 4') 

n'm'k' Pzi,Pz2,Pz3 
^^Pzc+Pzl,Pz2+Pz3^Pz,Pz\^nn' + ^p z ,Pz2^nm' ~ ^ P z,Pz3^nl') 



X 



(i + /r')/ 2 m 73 fc '-/r(i + / 2 m ')(i + /: 



(A13) 



Similarly, we can obtain the expression for C22 collision term, which is the contribution from 
the H' 4 perturbation. 

C22L/W] = -iTrp(t,to)[/(Pz,^to),^(t)] 



24 



—I 



y E E E ^ z 5(p zl + Pz2 ,p z3 +p z4 ) 

n 1 m 'le'l' 



rim'k'V qz Pzi,Pz2,Pz3 



X [dpzi,Pz+qz/2$rnn' ( a l z - gz /2n a p z 2m' a Pz3k' a p z 4l')t 
+^Pz2,Pz+qz/2^mm' - ?2 /2n a p z m' a p z 3 fe' °W ) * 



/ 2 8 n k> (a pzin > a Pz2 m> a Pz^' a p z +qz/2m)t 



J Pz3,Pz~qz 



$Pz4,Pz-qz/25nl> (a Pzl n' a p z 2m' a Pz3k'a> Pz +q z /2rn)t\ ■ 



(A14) 



Using Wick's theorem, we obtain the relevant contribution 



( a pzin a Pz 2m a Pz3k a P z4l)t — 27TZ(yf nm fc/(5(e 1 + 6 2 e 3 e A)^Pzi+Pz2,Pz3+Pz4 

x[/r/ 2 m (i + ; 3 fc )(i + /I) - (i + A n )(i + / 2 m )/ 3 Vi] 



(A15) 



Inserting (IA15[) into (IAI4p . we obtain 



C22 [/« 



7T 



E 

n'm'k'l 



9n 



'm'k'V E 

Pzl,Pz2,Pz3,Pz4 



6{ef 



~k' _ ~m' _ ~V- 
e 2 e 3 e 4, 



^-^Pzl+Pz2,Pz3+Pz4(^Pz,Pzl^nn' + ^p z ,Pz2^nk' ^Pz,Pz3^nm' ^Pz,Pz4^nl' , 



X 



/r7!'(i + / 3 m ')(i + Si) - (i + /i n ')(i + ti')f? A 



(A16) 



Appendix B: Numerical methods 

In this section we discuss solution of the collisionless Boltzmann equation using N-body 
simulations. The effect of collisions is dealt with later. It is generally very difficult to solve 
using standard methods for treating partial differential equations. An alternative approach 
used extensively in the literature is to represent the phase-space density f n (p z ,z,t) by a 
cloud of discrete test particles. The momentum and position of each particle in an external 
potential U n (z, t) is then evolved according to Newton's equations. The ith test particle has 
variables {zi(t) , pi(t) , rii(t)}. Test particles keep motion in one dimension along z-axis. 

The phase-space distribution for this situation is given by 

f n {Pz,z,t) = -~— 2vrft^5(> - Zi)S(p z -p zi )S n>ni . (Bl) 

where the weighting factor is fixed by the requirement that the phase-space distribution 
is normalized to the number of physical atoms, N, with NTi = J2nl dzdp z f n . By using a 
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sufficiently large number of test particles, Nt, a reasonable approximation to the continuous 
phase-space distribution is obtained. Note that the number of test and physical particles is 
not necessarily equal. In fact, for a relatively small number of physical atoms it is essential to 
simulate more test particles in order to minimize the effects of a discrete particle description. 

The time evolution of f n (p z ,z,t) is determined by the time-dependent potential and 
momentum variables of each test particle, given by the equations 



dzj _ p Zl {t) 
dt ~ M ' 

lould change by collision process. The collisions are treated in a similar manner 



The ni(t)s 

to Ref. p, ll9[ except calculation of angle. The ith test particle has the index of radial 
direction instead of angle. 

The phase-space variables are updated by advancing the position and momentum of each 
particle at discrete time steps At. Symplectic integrators are used extensively in molecular 
dynamics (MD) simulations since they possess several desirable properties, such as conser- 
vation of phase-space volume and of energy over a long period. We use a second-order 
symplectic integrator in our calculations, which is the classical analog of the split-operator 
method discussed earlier. To show this, it is convenient to work within the Lie formalism. 

2 

Consider the classical Hamiltonian for a single particle, Hi = + Vi(zi). The evolution of 
its phase-space coordinates Zi = (p zi , zi) is then determined by the equation 

^ = {Z i ,H i } = -i£Z i , (B3) 

where {F, G} = (d Zj Fd Pzi G — d Pzi Fd Zj Gj is the Poisson bracket and C is the Liouville 
operator. One can then write the time evolution Z as 

Z{t + At) = e- iCAt Z{t). (B4) 

Splitting the Hamiltonian into potential and kinetic terms, the effect of the classical operator 
in the simulations is to update the particle positions and velocities in three steps 

Zi = Zi + ~AtVi(t), 

v i (t + At)=v i -——U ni (z i ), 
M azi 

Zi(t + At) = Zi + -AtVi(t + At). (B5) 
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The effective potential U is determined self-consistently as the system evolves in time, 
and includes the condensate mean-field and the mean-field generated by the thermal cloud. 
The latter is in general much weaker than the condensate mean-field due to the larger 
spatial extent of the thermal cloud. Nevertheless, it is important to include this term in 
order to ensure the conservation of the total energy of the system. Although the calculation 
of the condensate mean-field is straightforward, the use of discrete particles with a contact 
interatomic potential creates a problem in determining the noncondensate mean-field. Taken 
literally, the mean-field consists of a series of delta peaks 

U2(z,t) = ^-2j2g nnkl J25(z-z i )5 n , ni = 2g n n^(z,t). (B6) 

J^tp k,l i=0 

This expression clearly cannot be used as it is to generate the forces acting on the test 
particles that are required in the MD simulation. We generate a smooth thermal cloud 
density by performing a convolution with a sampling (or smoothening) function S(z) which 
is normalized to unity. In particular, we define 

r ~ N 

U$( Z ,t) = / dz'S(z-z')U2( Z ',t) = ^Y.gnnklY.^-^nv (B7) 

where we choose S(z) ~ e~ z2 ^ 2 , i.e., an isotropic Gaussian sampling function of width 77. 

We proceed by making use of a FFT. First, each particle in the ensemble is assigned to 
points on the ID Cartesian grid using a cloud-in-cell method. We now consider a particle at 
position z, between two grid points at z k and z k+ ±. The particle is assigned to both points 
with weightings (1 — a) and a, respectively, where a = (z — z k )/(z k+1 — z k ). This can be 
viewed as a more sophisticated binning procedure in that it takes into account the actual 
positions of particles within the cells. We then convolve the cloud-in-cell density with the 
sampling function by Fourier transforming it and then multiplying it by the analytic FT of 
the sampling function. An inverse FFT then generates the sampled potential. This potential 
is used directly in the GP evolution, while the forces on the test particles are obtained by 
taking a numerical derivative and interpolating to the positions of the particles. We have 
also checked that small variations of rj about the value chosen to do the simulations have 
little effect on our final results. 

Probabilities for either C22 or C12 collisions are calculated in a way which is consistent with 
a Monte Carlo sampling of the collision integrals, as discussed below. We first give details 
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for the C22 integral, which physically corresponds to scattering of two thermal particles 
into two final thermal states. Hence the process conserves the number of thermal atoms 
/ dp z / (27ih)C 2 2 = 0. We are interested in the mean collision rate at a point z, which is given 
by 



r£U] = / ^2°2 U U], (B8) 

where 

C 2 ° 2 ut [/n] = E ^frr- I d Pz2 dp z3 dp z , - ef - ef ~ 

n'm'k'l' 47rft J 

x5{Pz + Pz2 - Pz3 - PzA){Snn' + S nk ' ~ S nm > - 5 nV ) 

x/r7 2 fe '(i+/3 m ')(i+/r)- (B9) 

We now write the required local collision rate as 

r^U] = E / It / If r'Mf*(pa)gS(p.i,P«) 

= E / rfp^ n ' fe '(p,)^(p,), (BIO) 

where p z is a point in two-dimensional momentum space and the factor 

w n ' k '( Pz ) = r\p zl )f k \p z2 )/{2nh)\ (Bll) 



is considered as a weight function. We denote the maximum value of w n (p z ) by w^!^ 
and define the domain on which the integrand is nonzero by [— p zlI1 ax/2,p zma , x /2] for each 
momentum component. Choosing a point p Z i at random in the hypervolume (p zmax ) 2 , and 
a random number Ri uniformly distributed on [0,w™'^] the point p Z i is is accepted if Ri < 
w^y- and the quantity g™2'<(Pzi) is accumulated. The value of the integral is then given 
approximately as 

n'm'k'l' iV i 

where N is the number of random p Z i points chosen and the prime on the summation includes 
only those points for which Ri < w nk (p iz ). For g%£, the integral is simply 

j\rnk 

h n (z)h k (z) = (p 2max ) 2 <^^, (B13) 



where N™ k is the total number of points accepted, and 
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rSf[/n]= E n n '(z)h k '(z)^j:'g^:(p zl ). (B14) 

n'm'k'l' 



The sample of iV™ fc points accepted consists of N™ k p zl values and N™ k p z2 values, each of 
which is distributed according to f n (p z \) and f k (p z2 ). This set of 2Ns nk p values can be 
identified with jV ce ii test particles in a cell of volume Az. If this set is to be representative 
of the local density, we must have 

n k (z) = — = (B15) 

With this identification, 

A*rsU]=2 £ ^'(z^gSipUM- ( B1 6) 

n'm'k'l' i 

In other words, the collision rate can be estimated by sampling the test particles in the cell 
Az in paris. 

For our purposes it is convenient to express the integral in terms of new momentum 
variables (p z0 ,p' z0 ) and (p" z ,p" z ): p z i, z2 = (p z o ± P")/VZ and p z3M = (p' z0 ±p z )/V2. p z0 
and p" are proportional to the center-of-mass and relative momenta, respectively, of the 
incoming 1 and 2 particles. The momentum and energy delta functions reduce to 

m + l\ - if - e{)5(p zl + Pz2 - p z3 - Pz4 ) = ^=5(p z0 - p' z0 )5(pf - pf + Kl)(B17) 
with E™1 = 2M(E n + E k — E t — E m ). Integrating over p' z0 and p' z , we obtain 

Tzft/n] = ^2 ~ 7 \ 2 (^tm' + &nk' ~ ^nm' ~ <W) 

n'm'k'l' 47r/1, 

/ dp zl ff J d p z2 f k '\p zl -p z2 \(l + / 3 m ')(l + fi). (B18) 

where p z3z 4 = P z o ± VP" 2 + ^nk Calculation of the rate therefore involves integrals over all 
possible initial states and all radial states. Inserting the explicit form of g for Y 22 collision 
rate, we have 

AzT° 2 f[f n ] = J2^Mg 2 nmkl J2n k (z)\p zi - p zj \(l + / 3 m )(l + /]), (B19) 
ki (ij) 

where the sum is now taken over pairs of test particles. This expression allows us to define 
the probability P?- 2 that a pair of atoms (ij) in the cell suffers a collision in a time interval 
At, 

Pfj'ifn] = nM 9l'm'k'l'($nn' + <W ~ Km' ~ Kl') 

fi k \z)\p zl -p Z3 \(l + fl n ')(l + fl). (B20) 
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Selecting atoms in pairs from each cell and assigning them a collision probability Pf 2 allows 
us to simulate the effect of collisions in a way which is consistent with the Boltzmann collision 
integral. 

We treat C\ 2 collisions somewhat differently. First, we note that the total rate of change 
of the number of thermal atoms per unit volume due to these collisions is 



/ ^k Cu[fn] = S g ^h° nc I dPz2 1 dPzs I dPz ^ ec + ^ - e ~™' - ~ e * )6{Pcz + Pz2 - Pz3 - PzA 



n'm'k 1 



X((W - Km' ~ <W) [^'(l + / 3 m ')(l + ft) - (1 + tf)ftfi 



(B21) 



According to this definition, using same transformation as in (IB19I) 

r?f [f n ] = £ ^mW0"c J ^ J ^ J dpz4§{€c + ~ e n_ ~m _ ~k } 
n'm'k' 

x (5 nn > - 5 nm , - 5 nk ,)5{j) cz + p z2 - p z3 - Pz^fZ'i 1 + •/•T'X 1 + /f) 
= E / d P J 2n ' m ' k '^ nc f? P r(l + / 3 m ')(l + f!')(5 nn , - 6 nm , - 6 nk ,), (B22) 

n'm'k' J 1X11 



where p° ut = y\p cz — Pz2\ 2 — ^M{U n + E n — /x c ). This rate can be estimated by writing 

r°u[fn] = E / dp Z 2W n '(p z2 )g^' kl ( Pz2 )(6 nnl - 5 nm , - 5 nk ,), (B23) 

n'm'k' 

where w n (p z2 ) = f n {p Z 2) I '{%kK) and g™ k {jpz2) is the remaining part of the integrand. A Monte 
Carlo sampling of the integral leads to the estimate 

N s 

AzT°£[Q ~ J! X>£(P*)(<W ~ Km' - <W), (B24) 

n'm'k' i=l 

where N s represents the number of atoms in the cell of volume Az. The probability of an 
atom in the cell suffering this kind of collision in the time interval At is therefore 



PZL[fn]=9nk(V Zl )At 

= E 2nC9 h k ° M ^\Pc. - PA 2 - + E n ~ He) 

nmk h 

x(l + / 3 m )(l + / 4 fc )At 

In r/ 2 M i 

= E lT ^\p cz -p Z 2\ 2 -4M(E n -E Q + g> nn00 n c ) 
k n 

x(l + / 3 m )(l + / 4 fe )At, (B25) 
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where we assume 

U n -fMc = V cxt (z) + 2(g nn00 n c + ^g nnkk h k ) - (V ext (z) + g 00 o n c + J2dookkn k ) 

k k 

— ^2( 2 9nnkk ~ gookk)n c = g' nn0Q n c . (B26) 
k 

The "in" collision rate is given by 

2 

= £ 9n ' m ' k J nc I d Pz2 [ dp z3 [ dp zi S(e c + 4 - ef - ~4){5 nn , - 6 nm , - 5 nk ,) 

n'm'k' Tlh J J J 

xS( Pcz + Pz2 - Pz3 - Pz a){1 + ft'Wfi 

Ef dp Z 2 fm > f dp Z 4 k i g n , m , k , Mn c fn'\/x x x \ 
^ J M h J M h — vn — (1 + h ){6nn ' - 6nm ' - 5nk,) 

n m k 

&{{Pcz - Pz±){Pcz - Pz2) - Mg' n , n , QO n c + M(E Q + E n > - E m , - E k >)], 

(B27) 

where we have interchanged the particle labels 2 and 3 to obtain the second line in this 
equation. This rate corresponds to two thermal atoms scattering into a condensate atom 
and an outgoing thermal atom, and is thus the rate that atoms feed into the condensate as 
a result of collisions. Although the collision of atoms 2 and 4 can be treated by the methods 
used to analyze the C22 collision rate, it is preferable to define a single atom collision rate 
by writing this integral in the form of Eq. ( 1B25I) and performing a Monte Carlo sampling 
with respect to the p z2 variable. This procedure leads to the collision probability per atom 
pin \f 1 _ a + f dp Z A fk g nmk0 Mn c . 

n nmk l/nj -M J —h — (1 + / s ) 

S[(pcz - Pz4)(Pcz - pU) - Mg' nn00 n c + M(E + E n - E m - E k )}. 

(B28) 

This analysis yields probabilities for a particular atom to undergo 'out' or 'in' collisions. To 
decide whether either event takes place, another random number < X 12 < 1 is chosen. 
If Xyi < P° nt then an 'out' collision is accepted; the incoming thermal atom is removed 
from the ensemble of test particles and two new thermal atoms are created. However, if 
pout < ^ 2 < pout _|_ pm ^ then an 'in' collision takes place and atom 2 is removed from the 
thermal sample. In addition, a second test particle, atom 4, is removed and a new thermal 
atom, atom 3, is created. In practice, it is exceedingly unlikely that a test particle will exist 
that will precisely match the required phase-space coordinates of particle 4. We therefore 
search for a test particle in neighboring phase-space cells and remove this particle if one is 
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found. This can be justified by remembering that we are only interested in describing the 
evolution in phase-space in a statistical wayis misleading to think of a direct correspondence 
between the test particles and physical atoms. If no test particle exists in the vicinity of v±, 
the local phase-space density f&, and hence P™, will be zero and the 'in' collision is precluded 
from occurring in any case. The above procedure leads to a change in the number of atoms 
in the thermal cloud. In order to conserve the total particle number the GP equation is 
propagated with the R term which changes the normalization of the wave function and hence 
the condensate number. This quantity can be evaluated from the Monte Carlo process decribed 
above by summing probabilities for particles 

= E Em?! - tJ- (B29) 

c nmk i 

In practice, this assignment to grid points is performed with a cloud-in-cell approach similar 
to the one described earlier. Of course, the normalization of the condensate wave function 
varies continuously as opposed to the variation of the thermal atom number which changes 
by discrete jumps. Nevertheless, one can show that the subsequent change in the condensate 
normalization is consistent with the addition or removal of atoms from the thermal cloud, 
so that the total particle number, N tot , is conserved within statistical fluctuations. 
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